#ifndef _BLES2_H_
#define _BLES2_H_

//#include "Equations/TwoPhaseEquation.hpp"
#include "./TwoPhaseEquation.hpp"

namespace Tuna {

  template <typename T, int Dim> class BLES2;

  template <typename T>
  struct Typeinfo <BLES2<T, 1> > {
    typedef T prec_t;
    enum { Dim = 1 };
  };

  template <typename T>
  struct Typeinfo <BLES2<T, 2> > {
    typedef T prec_t;
    enum { Dim = 2 };
  };

  template <typename T>
  struct Typeinfo <BLES2<T, 3> > {
    typedef T prec_t;
    enum { Dim = 3 };
  };
  
  template<typename Tprec, int Dim>
  class BLES2 : public TwoPhaseEquation<BLES2<Tprec, Dim> >
  {

    typedef TwoPhaseEquation<BLES2<Tprec, Dim> > TP_BLES2;

    using GeneralEquation< TP_BLES2 >::aE;
    using GeneralEquation< TP_BLES2 >::aW;
    using GeneralEquation< TP_BLES2 >::aN;
    using GeneralEquation< TP_BLES2 >::aS;
    using GeneralEquation< TP_BLES2 >::aF;
    using GeneralEquation< TP_BLES2 >::aB;
    using GeneralEquation< TP_BLES2 >::aP;
    using GeneralEquation< TP_BLES2 >::sp;
    using GeneralEquation< TP_BLES2 >::dx;
    using GeneralEquation< TP_BLES2 >::dy;
    using GeneralEquation< TP_BLES2 >::dz;
    using GeneralEquation< TP_BLES2 >::dt;
    using GeneralEquation< TP_BLES2 >::bi;
    using GeneralEquation< TP_BLES2 >::ei;
    using GeneralEquation< TP_BLES2 >::bj;
    using GeneralEquation< TP_BLES2 >::ej;
    using GeneralEquation< TP_BLES2 >::bk;
    using GeneralEquation< TP_BLES2 >::ek;
    using GeneralEquation< TP_BLES2 >::applyBoundaryConditions1D;
    using GeneralEquation< TP_BLES2 >::applyBoundaryConditions2D;
    using GeneralEquation< TP_BLES2 >::applyBoundaryConditions3D;

    using TP_BLES2::S;
    using TP_BLES2::p;
    using TP_BLES2::phi_0;
    using TP_BLES2::Srw;
    using TP_BLES2::Sro;
    using TP_BLES2::mu_w;
    using TP_BLES2::mu_o;
    using TP_BLES2::k;
    using TP_BLES2::porosity;

  public:
    typedef Tprec prec_t;
    typedef typename TunaArray<prec_t, Dim >::huge ScalarField;
    
    BLES2() : TwoPhaseEquation<BLES2<prec_t, Dim > >() { }    
    ~BLES2() { };
    
    inline bool calcCoefficients1D(); 
    inline bool calcCoefficients2D();
    inline bool calcCoefficients3D() { };
    inline void printInfo() { std::cout << " BLES2 "; }
  };

/*
 *  Quadratic for realtive permeability , Upwind for Sw
 */
template<typename Tprec, int Dim>
inline bool BLES2<Tprec, Dim>::calcCoefficients1D () 
{
    static prec_t Sw_e, Sw_w;

    // Quadratic
    static prec_t mult = k * dt / 
      ( porosity * dx * dx * (1 - Srw - Sro) * (1 - Srw - Sro) * mu_w );
      
    aE = 0.0; aW = 0.0; aP = 0.0; sp = 0.0;

    for (int i =  bi; i <= ei; ++i) {      

      // Upwind
      if ( phi_0(i+1) > phi_0(i) ) Sw_e = phi_0(i+1);
      else                          Sw_e = phi_0(i);
      if ( phi_0(i-1) > phi_0(i) ) Sw_w = phi_0(i-1);
      else                          Sw_w = phi_0(i);

      // Quadratic
      aE (i) = (Sw_e - Srw) * (Sw_e - Srw) * mult;
      aW (i) = (Sw_w - Srw) * (Sw_w - Srw) * mult;
      aP (i) = aE (i) + aW (i);      
    }

/* ----- Boundary conditions ----- */
    // Dirichlet right side
    aP(ei) += aE(ei);
    sp(ei) = 2 * aE(ei) * p(ei+1);
    aE(ei) = 0;

    // Neumann left side
    aP(bi) -= aW(bi) ;
    sp(bi) = aW(bi) * dx * ( 3.47e-7 * mu_w / k) ;
    aW(bi) = 0;

    return 0;
}


/*
 *  Quadratic for realtive permeability , Upwind for Sw
 */
template<typename Tprec, int Dim>
inline bool BLES2<Tprec, Dim>::calcCoefficients2D () 
{
    static prec_t Sw_e, Sw_w, Sw_n, Sw_s;

    // Quadratic
    static prec_t multx = k * dt / 
      ( porosity * dx * dx * (1 - Srw - Sro) * (1 - Srw - Sro) * mu_w );

    static prec_t multy = k * dt / 
      ( porosity * dy * dy * (1 - Srw - Sro) * (1 - Srw - Sro) * mu_w );
      
    aE = 0.0; aW = 0.0; aN = 0.0; aS = 0.0; aP = 0.0; sp = 0.0;


    for (int i =  bi; i <= ei; ++i)    
      for (int j = bj; j <= ej; ++j) {  

	// Upwind
	if ( phi_0(i+1, j) >= phi_0(i, j) ) Sw_e = phi_0(i+1, j);
	else                                Sw_e = phi_0(i  , j);
	if ( phi_0(i-1, j) >= phi_0(i, j) ) Sw_w = phi_0(i-1, j);
	else                                Sw_w = phi_0(i  , j);
	
	if ( phi_0(i, j+1) >= phi_0(i, j) ) Sw_n = phi_0(i, j+1);
	else                                Sw_n = phi_0(i, j);
	if ( phi_0(i, j-1) >= phi_0(i, j) ) Sw_s = phi_0(i, j-1);
	else                                Sw_s = phi_0(i, j);

	// Quadratic
	aE (i, j) = (Sw_e - Srw) * (Sw_e - Srw) * multx;
	aW (i, j) = (Sw_w - Srw) * (Sw_w - Srw) * multx;
	aN (i, j) = (Sw_n - Srw) * (Sw_n - Srw) * multy;
	aS (i, j) = (Sw_s - Srw) * (Sw_s - Srw) * multy;	
	aP (i, j) = aE (i, j) + aW (i, j) + aN (i, j) + aS (i, j);
      }

  /* ----- Boundary conditions ----- */

    applyBoundaryConditions2D();
    
    // I'm using an explicit scheme, so the boundary conditions are
    // different than in the implicit case.
    
    Range I(bi,ei), J(bj, ej);
    
    // Dirichlet right side
    aP(ei, J) += aE(ei, J);
    sp(ei, J) = 2 * aE(ei, J) * p(ei+1, J);
    aE(ei, J) = 0;
    
    // Neumann left side
    // Remember that Sw_e = 0.8 on the boundary, that is
    // the reason why sp(bi) has the form below.
    aP(bi, J) -= aW(bi, J) ;
    sp(bi, J) = aW(bi, J) * dx * ( 3.47e-7 * mu_w / k) ;
    aW(bi, J) = 0;
    
    return 0;
}


} // Tuna namespace


#endif //_BLES2_H_

















